rm(list = ls())
cat("\014")

library(tidyverse, quietly = T)
library(sf,        quietly = T)
library(gganimate, quietly = T)
library(openair,   quietly = T)
library(lubridate, quietly = T)
library(zoo,       quietly = T)
library(ggthemes,  quietly = T)
library(ggmap,     quietly = T)
library(sp,        quietly = T)
library(magick,    quietly = T)

Import the time-activity diary for the line and tidy it up.

diary      <- read_csv('air_quality_data/victoria_line_2014-12-08.csv', col_types = cols())

base       <- tibble(date_time = seq(min(diary$start_time), max(diary$end_time), by = 'min'))

base       <- left_join(base, diary, by = c('date_time' = 'start_time')) %>% 
              select(-fake_id)

base       <- left_join(select(base, -end_time),
                        filter(base, date_time != end_time),
                        by = c('date_time' = 'end_time')) %>%
              mutate(line.x            = if_else(is.na(line.x)            & !is.na(line.y),            line.y,            line.x),
                     station.x         = if_else(is.na(station.x)         & !is.na(station.y),         station.y,         station.x),
                     environment.x     = if_else(is.na(environment.x)     & !is.na(environment.y),     environment.y,     environment.x),
                     sub_environment.x = if_else(is.na(sub_environment.x) & !is.na(sub_environment.y), sub_environment.y, sub_environment.x)) %>%
              select(date_time, line.x, station.x, environment.x, sub_environment.x) %>%
              rename(line = line.x, station = station.x, environment = environment.x, sub_environment = sub_environment.x)

rm(diary)

Import the measured air quality to join to the time-activity diary

air_quality <- read_csv('air_quality_data/LUexportMay2016.csv', col_types = cols()) %>%
                rename_all(tolower) %>%
                mutate(datetime = as.POSIXct(datetime, format='%d/%m/%Y %H:%M')) %>%
                filter(sitecode == 'CAR')

Join the time activity and air quality

line       <-  left_join(base, air_quality, by = c('date_time' = 'datetime')) %>% 
                mutate(scaledvalue = if_else(species == 'PM25', value, scaledvalue)) %>%
                select(-value) %>%
                rename(concentration = scaledvalue)

rm(air_quality, base)

Import the background air quality for that day to correct the PM2.5 measurements

background_pm25 <- importKCL(site = "kc1", year = c(2014,2015,2016), pollutant = "pm25", met = FALSE,
                             units = "mass", extra = FALSE)
## NOTE - mass units are used 
## ug/m3 for NOx, NO2, SO2, O3; mg/m3 for CO
## PM10_raw is raw data multiplied by 1.3
background_pm25 <- data.frame(background_pm25, day = as.Date(format(background_pm25$date)))
background_pm25 <- aggregate(pm25 ~ day, background_pm25, mean)

Join the background PM2.5 to the monitored concentrations

line     <-    mutate(line, day = date(date_time)) %>%
                left_join(background_pm25, by = c('day' = 'day')) %>%
                select(-day)

Correct the PM2.5 data

line <-         mutate(line, concentration = 
                         case_when(concentration >  pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ pm25 + (1.82 * (concentration - (pm25/0.44))),
                                   concentration <= pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ concentration * 0.44,
                                   TRUE ~ concentration)) %>%
                select(-pm25)

Get the station locations

stations <- st_read('https://raw.githubusercontent.com/dracos/underground-live-map/master/bin/stations.kml', quiet = T) %>%
            select(-Description) %>%
            rename_all(tolower) %>%
            mutate(name = gsub(' Station', '', name)) %>%
            mutate(name = gsub("'", '', name)) %>%
            st_zm(drop = TRUE)

Join locations to the concentrations and station names

line   <-left_join(line, stations, by = c('station' = 'name')) %>% st_as_sf()

Fill in the missing coordinates using linear interpolation

st_geometry(line) <- as_tibble(st_coordinates(line)) %>% 
                      mutate(X = na.approx(X), Y = na.approx(Y)) %>% 
                      st_as_sf(coords = c('X', 'Y')) %>%
                      st_geometry()

Start trying to animate. First the graph animation.

text_placement <- filter(line,species == 'PM25') %>%
                    as_tibble() %>%
                    select(-geometry) %>%
                    group_by(species) %>%
                    summarise(max_x = max(date_time), max_y = max(concentration))

line           <- left_join(line, text_placement, by = c('species' = 'species'))

pm25 <- filter(line,species == 'PM25' & !is.na(concentration))

pm25$moving_mean <- NA
pm25$pc25        <- NA
pm25$pc75        <- NA

for (i in 2:nrow(pm25)) { 
  pm25[i,'moving_mean'] <- as.integer(round(mean(pm25[1:i,]$concentration),0))
  pm25[i,'pc25']        <- as.integer(round(quantile(pm25[1:i,]$concentration, 0.25),0))
  pm25[i,'pc75']        <- as.integer(round(quantile(pm25[1:i,]$concentration, 0.75),0))
  }

graph <- ggplot(data = filter(line,species == 'PM25')) +
    geom_path(aes(date_time, concentration), colour = '#0099CC', size = 1.2) +
    geom_ribbon(data = filter(pm25,!is.na(pc25)), aes(x=date_time, ymin = pc25, ymax = pc75), alpha = 0.2) +
    geom_hline(data = filter(pm25,!is.na(pc25)), aes(yintercept = moving_mean), alpha = 0.2, colour = 'red') +
    geom_point(aes(date_time, concentration), colour = 'red', size = 4) +
    geom_text(data = filter(pm25,!is.na(pc25)), aes(x = max(date_time), y = moving_mean, label = moving_mean), hjust = 0, size = 8, colour = 'red') +
   geom_vline(xintercept = as.POSIXct('2014-12-08 10:00'), size = 1) +
    coord_cartesian(clip = 'off') +
    geom_text(aes(min(date_time), max_y,label = round(concentration,0)),
              size = 10, hjust = 0, vjust = 1, colour = '#0099CC') + 
    theme(axis.text.x  = element_blank(),
          axis.text.y  = element_text(size = 12),
          axis.title.y = element_text(size = 12),
          axis.line    = element_line(colour = 'black'),
          plot.margin  = margin(5.5, 40, 5.5, 5.5)) +
    xlab('') +
    ylab('PM2.5 ug/m3') +
    transition_reveal(date_time)

graph_animation <- animate(graph, fps = 3, height = 600, renderer = gifski_renderer(loop = FALSE))
anim_save(file='graph_animation.gif', animation = graph_animation, path = 'outputs')

graph_animation

Now the spatial animation.

line <- st_set_crs(line,4326)

map_area <- c(as.vector(st_bbox(line))[1]-0.01, as.vector(st_bbox(line))[2]-0.01, as.vector(st_bbox(line))[3]+0.01, as.vector(st_bbox(line))[4]+0.01)

map <-  ggmap(get_map(location = map_area), source = "google", maptype = "roadmap", zoom = 13) +
    geom_path(data = as_tibble(st_coordinates(filter(line,species == 'PM25')[1:56,])), aes(X, Y), colour = "#0099CC", size = 1) +
     geom_sf(data = filter(line,species == 'PM25'), size=4, colour = 'red',inherit.aes = FALSE) +
     coord_sf(datum=NA)  +
     transition_time(date_time) +
     ease_aes('linear') +
     theme(axis.title = element_blank(),
           panel.background = element_blank(),
           panel.border = element_rect(colour = 'black', fill = NA),
           axis.text = element_blank())

map_animation  <- animate(map, fps = 3, height = (470*2), width = (200*2))

anim_save(file='map_animation.gif', animation = map_animation, path = 'outputs')

map_animation